Basic Analyses of Isotopic Variation

In this practical, we will explore our data in more detail using some basic statistical tests to identify where the variation in our data lies, and hence where the interesting ecological insights might lie. In this example dataset, we have d13C, d15N and C:N ratio from three fish species.

mydata <- read.csv("Practical02.csv", header = TRUE, stringsAsFactors = FALSE)

# force species to be a factor type variable
mydata$Species <- factor(mydata$Species)

# verify that it imported as expected
head(mydata)

NA

Lipid correction

# code to do this using base R graphics
# plot(Liver.d13C ~ Liver.CN, data = mydata, 
#      subset = Species == "Arrowtooth.Flounder", main = "Arrowtooth.Flounder")
# 
# # lowess smoother not working like this for the subset.
# # lines(lowess(mydata$Liver.CN, mydata$Liver.d13C), col = "red")
# 
# plot(Liver.d13C ~ Liver.CN, data = mydata, 
#      subset = Species == "Pacific.cod", main = "Pacific.cod")
# 
# plot(Liver.d13C ~ Liver.CN, data = mydata, 
#      subset = Species == "Pollock", main = "Pollock")

# try to do this in ggplot using facets
library(tidyverse)
liver.plot <- ggplot(data = mydata, 
                     mapping = aes(x = Liver.CN, 
                                   y = Liver.d13C)) + 
  geom_point() +
  geom_smooth(method = "loess", alpha = 0.5) +
  facet_grid(Species~.) + 
  theme(text = element_text(size=20))
print(liver.plot)

`geom_smooth()` using formula 'y ~ x'

The liver in particular has a large quantity of lipids, and in this example, we apply a simple numerical correction based on Logan et al 2008. This correction is weighted by the C/N ratio, such that high amounts of C require a greater correction.

# constants used in the correction for liver tissue
a <-   6.059
b <- -22.270
c <-  -1.397

# add columns for lipid corrected liver to the data.frame
mydata <- mydata %>% 
  mutate(Liver.d13C.crtd = Liver.d13C + 
           (a * Liver.CN + b) / (Liver.CN + c))

# mydata$Liver.d13C.crtd <-  mydata$Liver.d13C + 
#   (a * mydata$Liver.CN + b) / (mydata$Liver.CN + c)

# check that it worked as expected
head(mydata)

# specify a color pallete with 3 colours for this example
palette(c("black", "red", "blue"))

# plot the correlation between the corrected and raw values, and 
# indicate which observations have large C/N ratios by the 
# size of the datapoints.
# plot(Liver.d13C.crtd ~ Liver.d13C, data = mydata, 
#      col = Species, cex = 5 * Liver.CN / max(Liver.CN), pch = 19, asp = 1)
# 
# # add the 1:1 line
# abline(a = 0, b = 1, col = "grey", lty = 1, lwd = 2)
# 
# 
# # reset the colour palette to the default
# palette("default")

# now plot the corrected values on top of our raw data from before 
# and add a vertical line at 3.6 which is "pure protein"
liver.plot.2 <- liver.plot + 
  geom_point(data = mydata, aes(x = Liver.CN, 
                                y = Liver.d13C.crtd), 
             col = "red") +
  geom_vline(xintercept = 3.6, col = "black")

print(liver.plot.2)

`geom_smooth()` using formula 'y ~ x'

Basic summary statistics

As always, it is sensible to start by taking a quick look at the summary statistics describing the data, and how the data are encoded in R. The only bit of code I’ve changed here is to make the text of the correlation coefficient text a bit smaller when defining the function panel.cor().

# how are the data encoded in R?
str(mydata)

'data.frame':   113 obs. of  9 variables:
 $ Species        : Factor w/ 3 levels "Arrowtooth.Flounder",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ Length         : int  320 170 180 360 160 500 420 350 370 480 ...
 $ Muscle.d13C    : num  -19.1 -19.3 -19.3 -18.5 -19.5 ...
 $ Muscle.d15N    : num  15.3 13.3 12.8 15.3 13 ...
 $ Muscle.CN      : num  3.21 3.21 3.13 3.23 3.27 ...
 $ Liver.d13C     : num  -22.5 -21.4 -22.2 -20.1 -22.6 ...
 $ Liver.d15N     : num  15.1 13 13.2 14.5 12.3 ...
 $ Liver.CN       : num  6.96 7.5 8.65 4.43 9.14 ...
 $ Liver.d13C.crtd: num  -18.9 -17.6 -18.1 -18.6 -18.4 ...

# a basic summary of each column in the data.frame
summary(mydata)

                Species       Length     Muscle.d13C      Muscle.d15N   
 Arrowtooth.Flounder:29   Min.   : 90   Min.   :-23.41   Min.   :12.49  
 Pacific.cod        :39   1st Qu.:270   1st Qu.:-19.75   1st Qu.:14.57  
 Pollock            :45   Median :480   Median :-18.83   Median :16.03  
                          Mean   :450   Mean   :-19.05   Mean   :16.18  
                          3rd Qu.:620   3rd Qu.:-18.11   3rd Qu.:17.79  
                          Max.   :820   Max.   :-16.77   Max.   :20.10  
   Muscle.CN        Liver.d13C       Liver.d15N        Liver.CN     
 Min.   : 3.105   Min.   :-28.48   Min.   : 9.896   Min.   : 3.225  
 1st Qu.: 3.184   1st Qu.:-24.08   1st Qu.:13.336   1st Qu.: 6.395  
 Median : 3.235   Median :-22.71   Median :15.125   Median :10.189  
 Mean   : 3.901   Mean   :-22.62   Mean   :15.032   Mean   :15.793  
 3rd Qu.: 3.327   3rd Qu.:-21.38   3rd Qu.:16.750   3rd Qu.:23.063  
 Max.   :21.149   Max.   :-17.35   Max.   :20.273   Max.   :70.376  
 Liver.d13C.crtd 
 Min.   :-22.85  
 1st Qu.:-19.03  
 Median :-18.41  
 Mean   :-18.62  
 3rd Qu.:-17.92  
 Max.   :-16.23  

# how many obverations do we have per Taxon might also be of use
table(mydata$Species)

Arrowtooth.Flounder         Pacific.cod             Pollock 
                 29                  39                  45 

# Using datq from our first practical, generate all pair-wise scatterplots
# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
# source("my_functions.R")
# These functions are copied from the help file for pairs()
## put histograms on the diagonal
panel.hist <- function(x, ...)
{
    usr <- par("usr"); on.exit(par(usr))
    par(usr = c(usr[1:2], 0, 1.5) )
    h <- hist(x, plot = FALSE)
    breaks <- h$breaks; nB <- length(breaks)
    y <- h$counts; y <- y/max(y)
    rect(breaks[-nB], 0, breaks[-1], y, col = "cyan", ...)
}

panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...)
{
    usr <- par("usr"); on.exit(par(usr))
    par(usr = c(0, 1, 0, 1))
    r <- cor(x, y)
    txt <- format(c(r, 0.123456789), digits = digits)[1]
    txt <- paste0(prefix, txt)
    if(missing(cex.cor)) cex.cor <- 0.4/strwidth(txt)
    # text(0.5, 0.5, txt, cex = cex.cor * r)
    text(0.5, 0.5, txt, cex = 2)
}

# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =

# this select function comes from the package dplyr and allows us to 
# select a subset of columns to plot.
pairs(mydata, 
      diag.panel = panel.hist,
      upper.panel = panel.smooth,
      lower.panel = panel.cor)

NA
NA
NA

Focus on variation among species

Using code from the first practical, we can use the package dplyr to generate summary statistics on the means and variation of each species in the dataset. I summarise the muscle and liver (using the lipid-corrected d13C) data separately only to keep the resulting table smaller and nicer to print.

# Summarise the data by group using the package dplyr (install if required)

# load the library
library(dplyr)

summarise.muscle <- mydata %>% group_by(Species) %>%  
  summarise(count = length(Species),
            musC = mean(Muscle.d13C), musSdC =  sd(Muscle.d13C), 
            musN = mean(Muscle.d15N), musSdN =  sd(Muscle.d15N),
            musCN = mean(Muscle.CN),  musSdCN = sd(Muscle.CN))

# print a pretty table
knitr::kable(summarise.muscle, digits = 2)

Species	count	musC	musSdC	musN	musSdN	musCN	musSdCN
Arrowtooth.Flounder	29	-20.08	1.30	14.64	1.2	4.30	2.21
Pacific.cod	39	-18.07	1.07	17.98	1.0	3.78	2.41
Pollock	45	-19.23	1.12	15.60	1.7	3.75	2.82

# print(summarise.muscle)


summarise.liver <- mydata %>% group_by(Species) %>% 
  summarise(count = length(Species),
            livC.c = mean(Liver.d13C.crtd), 
            livSdC.c =  sd(Liver.d13C.crtd),
            livN = mean(Liver.d15N), 
            livSdN =  sd(Liver.d15N),
            livCN = mean(Liver.CN),  
            livSdCN = sd(Liver.CN) )

# print a pretty table
knitr::kable(summarise.liver, digits = 2)

Species	count	livC.c	livSdC.c	livN	livSdN	livCN	livSdCN
Arrowtooth.Flounder	29	-18.66	0.70	14.21	1.34	8.38	3.79
Pacific.cod	39	-18.02	1.18	16.77	1.48	12.65	7.97
Pollock	45	-19.12	1.22	14.06	1.97	23.29	15.43

# print(summarise.liver)

The most suitable way to visualise the variation in the key tissues by species is probably to use simple boxplots. From here on in we will probably focus on the muscle data only and will put analysis of the liver samples on hold.

boxplot(Length ~ Species, data = mydata, col="grey", 
        xlab="Species", ylab="Length (cm)", 
        main = " variation in fish size")

boxplot(Muscle.d13C ~ Species, data = mydata, col="grey", 
        xlab="Species", ylab="Muscle d13C", 
        main = " variation in muscle d13C")

#now d15N
boxplot(Muscle.d15N ~ Species, data = mydata, col="grey", 
        xlab="Species", ylab="Muscle d15N", 
        main = " variation in muscle d15N")

#now d15N
boxplot(log(Muscle.CN) ~ Species, data = mydata, col="grey", 
        xlab="Species", ylab="log Muscle C/N", 
        main = " variation in muscle C/N")

NA
NA

An appropriate statistical test to determine whether the means of these boxplots are significantly different from one another is a simple ANOVA. You need to take a little care here, as the function you want when analysing data is called aov() in R (for Analysis of Variance) whereas there is also a function anova() but it is used to compare fitted linear model objects (a discussion point for another day, and one i wouldnt recommend in any case). I have log-transformed the C/N data owing to the heavy skew evident in the pairs plot above.

# length by species
length.aov <-  aov(Length ~ Species, data = mydata)
summary(length.aov)

             Df  Sum Sq Mean Sq F value   Pr(>F)    
Species       2  572074  286037   7.872 0.000638 ***
Residuals   110 3996926   36336                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# muscle carbon by species
muscleC.aov <- aov(Muscle.d13C ~ Species, data = mydata)
summary(muscleC.aov)

             Df Sum Sq Mean Sq F value   Pr(>F)    
Species       2  69.47   34.74   26.07 5.41e-10 ***
Residuals   110 146.57    1.33                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# muscle N by species
muscleN.aov <- aov(Muscle.d15N ~ Species, data = mydata)
summary(muscleN.aov)

             Df Sum Sq Mean Sq F value Pr(>F)    
Species       2  210.8  105.38   56.49 <2e-16 ***
Residuals   110  205.2    1.87                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# muscle C/N ratio by species
muscleCN.aov <- aov(log(Muscle.CN) ~ Species, data = mydata)
summary(muscleCN.aov)

             Df Sum Sq Mean Sq F value Pr(>F)
Species       2  0.437  0.2184   2.044  0.134
Residuals   110 11.757  0.1069               

Identifying which species are different from each other can be acheived using post-hoc tests that correct the p-values to take account of the multiple testing involved. Just looking at the d13C muscle data in this example (you can try the others as you like):

TukeyHSD(muscleC.aov)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = Muscle.d13C ~ Species, data = mydata)

$Species
                                      diff        lwr        upr     p adj
Pacific.cod-Arrowtooth.Flounder  2.0063584  1.3338843  2.6788325 0.0000000
Pollock-Arrowtooth.Flounder      0.8457546  0.1926797  1.4988296 0.0073897
Pollock-Pacific.cod             -1.1606038 -1.7606068 -0.5606007 0.0000343

Is it length or species that drives variation in d13C and d15N?

The larger Pacific Cod appears to be the largest of the fish on average, and they also have the larget d13C and d15N values. We might rightly want to know if this is due to their larger body size (length) or whether it is true species variation. We use analysis of covariance to investigate.

# specify a color pallete with 3 colours for this example
palette(c("black", "red", "blue"))

par(mfrow = c(1,2))
plot(Muscle.d13C ~ Length, col = mydata$Species, data = mydata, 
     main = "Muscle d13C")
legend("topleft", levels(mydata$Species), col = 1:3, pch = 1)

plot(Muscle.d15N~ Length, col = mydata$Species, data = mydata,
     main = "Muscle d15N")
legend("topleft", levels(mydata$Species), col = 1:3, pch = 1)

# reset the palette
palette("default")

d13C.length.model <- glm(Muscle.d13C ~ Species + Length, data = mydata)
summary(d13C.length.model)

Call:
glm(formula = Muscle.d13C ~ Species + Length, data = mydata)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.8335  -0.3881   0.3021   0.7082   1.7101  

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)        -2.024e+01  2.867e-01 -70.593  < 2e-16 ***
SpeciesPacific.cod  1.927e+00  2.981e-01   6.464 2.96e-09 ***
SpeciesPollock      7.635e-01  2.913e-01   2.621    0.010 *  
Length              4.975e-04  5.781e-04   0.861    0.391    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 1.335639)

    Null deviance: 216.05  on 112  degrees of freedom
Residual deviance: 145.58  on 109  degrees of freedom
AIC: 359.31

Number of Fisher Scoring iterations: 2

d15N.length.model <- glm(Muscle.d15N ~ Species + Length, data = mydata)
summary(d15N.length.model)

Call:
glm(formula = Muscle.d15N ~ Species + Length, data = mydata)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.7003  -1.0161   0.1646   0.8171   3.4259  

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)        1.375e+01  3.152e-01  43.629  < 2e-16 ***
SpeciesPacific.cod 2.910e+00  3.277e-01   8.880 1.52e-14 ***
SpeciesPollock     5.111e-01  3.202e-01   1.596    0.113    
Length             2.707e-03  6.354e-04   4.260 4.35e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 1.613847)

    Null deviance: 415.96  on 112  degrees of freedom
Residual deviance: 175.91  on 109  degrees of freedom
AIC: 380.69

Number of Fisher Scoring iterations: 2

While d15N appears to scale positively with length (and you can probably think why this might be), d13C does not appear to be affected by length. You could assess this by comparing a model without length using AIC, where we are looking for AIC to be lower by more than 2 units if we are to justify the inclusion of length in the model:

d13C.species.model <- glm(Muscle.d13C ~ Species, data = mydata)
summary(d13C.species.model)

Call:
glm(formula = Muscle.d13C ~ Species, data = mydata)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.6913  -0.4828   0.2836   0.7510   1.6583  

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)        -20.0782     0.2144 -93.668  < 2e-16 ***
SpeciesPacific.cod   2.0064     0.2830   7.088 1.37e-10 ***
SpeciesPollock       0.8458     0.2749   3.077  0.00264 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 1.332492)

    Null deviance: 216.05  on 112  degrees of freedom
Residual deviance: 146.57  on 110  degrees of freedom
AIC: 358.08

Number of Fisher Scoring iterations: 2

Based on this comparison, with an AIC of 359.31 for the model with length, and an AIC of 358.08 without length, we can conclude that length does not add sufficient explanatory information to warrant its inclusion.

Task: You might want to check whether there is an interaction between species and length in this relationship. Essentially this involves fitting 3 lines, one for each species by length. You can add an interaction to the glm() formulae, using a : to specify which variables interact.

Pair-wise comparisons of the different tissues

We might well ask whether the isotope values are the same between the liver and muscle tissues within the individual fish. To acheive this, we can conduct a pair-wise t-test that adds power over a standard t-test by treating the observations as being matched, or paired. Since we can’t use the formula method (y~x) easily when we have two columns (vectors) of data to compare, we can’t use the data = mydata option, and so we have to use the $ notation to provide the two vectors for comparison. We can’t even use the subset = option so instead we have to subset each of our vectors manually: mydata$Muscle.d13C[mydata$Species == "Arrowtooth.Flounder"].

You could take this code, copy it and modify it to run paired t-tests on the other species.

t.test(mydata$Muscle.d13C[mydata$Species == "Arrowtooth.Flounder"],
       mydata$Liver.d13C.crtd[mydata$Species == "Arrowtooth.Flounder"],
       paired = TRUE)

    Paired t-test

data:  mydata$Muscle.d13C[mydata$Species == "Arrowtooth.Flounder"] and mydata$Liver.d13C.crtd[mydata$Species == "Arrowtooth.Flounder"]
t = -6.1193, df = 28, p-value = 1.328e-06
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -1.8971660 -0.9455698
sample estimates:
mean of the differences 
              -1.421368 

Use multivariate statistics to compare d13C and d15N simultaneously among species

Visualising the multivariate data is easiest if we add ellipses over each species group as we did in the first practical session. The code that follows is taken from there and modified accordingly to make sure the \(x\) and \(y\) data names are appropriate for this new dataset.

Now we can make a nice plot with ggplot taking care to change the aesthetics to make our column names if we just copied and pasted this in from the first practical file.

# load the library
library(ggplot2)

# use our ellipse function to generate the ellipses for plotting

# this is the basic plot, which im not actually going to plot
# as i won't call print(first.plot). Doing this is purely choice, 
# and for me it keeps my code tidier and easier to interpret
first.plot <- ggplot(data = mydata, aes(Muscle.d13C, Muscle.d15N)) +
  geom_point(aes(color = Species), size = 2)+
  ylab(expression(paste(delta^{15}, "N (\u2030)")))+
  xlab(expression(paste(delta^{13}, "C (\u2030)"))) + 
  theme(text = element_text(size=15))

# decide how big an ellipse you want to draw
# NB 50% ellipses this time for no reason other than i dont need them huge
# to get a sense for their size and shape.. indeed we could plot 
# Standard Ellipses using p.ell <- stats::pchisq(1, df = 2) which
# results in 0.39
p.ell <- 0.50


# create our plot based on first.plot above
# adding the stat_ellipse() geometry. We 
# specify thee ellipse to be plotted using 
# the polygon geom, with fill and edge colour
# defined by Taxon as a grouping variable, 
# using the normal distribution and with 
# a quite high level of transparency.
ellipse.plot <- first.plot + 
  stat_ellipse(aes(group = Species, 
                   fill = Species, 
                   color = Species), 
               alpha = 0.3, 
               level = p.ell,
               type = "norm",
               geom = "polygon")

print(ellipse.plot)

Analysing these data to test whether the means of both the d13C and the d15N muscle data are simultaneously different can be achieved using Multivariate Analysis of Variance (MANOVA) or PERMANOVA if you are not happy with the assumptions of MANOVA, e.g. that the data are multivariate normal and that they have similar (assumed the same) covariance. In this case, both assumptions seem pretty reasonable to me so I will go with MANOVA here in the first instance. Both models require that we bind the two (or more) column vectors of data together into a matrix which we can do using cbind(). Remember that both these methods take the subet = option should you wish to restrict your comparison to a subset of groups. A useful trick here might be to use the match function (which has an odd, but nice to use alias %in%) and in this example something like: subset = (Species %in% c("Arrowtooth.Flounder", "Pacific.cod")) could provide a useful template for you, especially should you wish to restrict a larger dataset to only the fish or only the invertebrates if they were all in the same dataset. Alternateively, you could create new data.frames of only the data you want, and call them something helpful; you could acheive this with the function dplyr::filter()

multivar.model <- manova(cbind(Muscle.d13C, Muscle.d15N) ~ Species, 
                         data = mydata)

summary(multivar.model)

           Df  Pillai approx F num Df den Df    Pr(>F)    
Species     2 0.53826   20.253      4    220 3.223e-14 ***
Residuals 110                                             
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Permanova

Running the PERMANOVA analysis requires the package vegan so you will need to install.packages("vegan") if you do not already have it. The PERMANOVA method is a randomisation process that jumbles up the data among the groups and determines how likely our observed data are given random chance. It therefore does make the same parametric assumptions about the data being multivariate normal distributed with a common covariance structure among groups.

library(vegan)

Loading required package: permute
Loading required package: lattice
This is vegan 2.5-7

# extract the isotope data for muscle for 
# subsequent modelling as a response variable.
Y_muscle <- with(mydata, cbind(Muscle.d13C, Muscle.d15N))

# run a PERMANOVA model
perm.model <- adonis(Y_muscle ~ Species, 
                     data = mydata, 
                     method = "euclidean", 
                     permutations = 9999)

# print output of the permanova model to screen
perm.model

Call:
adonis(formula = Y_muscle ~ Species, data = mydata, permutations = 9999,      method = "euclidean") 

Permutation: free
Number of permutations: 9999

Terms added sequentially (first to last)

           Df SumsOfSqs MeanSqs F.Model     R2 Pr(>F)    
Species     2    280.23 140.117  43.815 0.4434  1e-04 ***
Residuals 110    351.77   3.198         0.5566           
Total     112    632.01                 1.0000           
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Distance between centroids

Some additional context we might like to add to this analysis is how far apart the centroids are. The centroids of these ellipses are simply the mean(d13C) and the mean(d15N) of each group, and we have already calculated these for the graph above and the results held in the object summarise.muscle. We can then extract the two columns of means, and use the dist() function in R to calculate the pairwise euclidean distances between them.

centroids <- cbind(summarise.muscle$musC, summarise.muscle$musN)

dist.btw.centoids <- dist(centroids)

print(dist.btw.centoids)

         1        2
2 3.898797         
3 1.278276 2.651899

Use the betadisper method to peform pairwise comparisons

# calculate the euclidean distances between observations
dis <- vegdist(Y_muscle, 
               method ="euclidean")

# analysis of multivariate homogeneity of 
# group dispersions (variances)
mod <- betadisper(dis, mydata$Species)

# print the results of the ANOVA to screen
knitr::kable(anova(mod))

	Df	Sum Sq	Mean Sq	F value	Pr(>F)
Groups	2	7.336073	3.6680365	4.79097	0.0101156
Residuals	110	84.217605	0.7656146	NA	NA

# perform Tukey's posthoc test between the means
mod.HSD <- TukeyHSD(mod)

# plot the post hoc test results
par(mfrow = c(1,1))
plot(mod.HSD)

Perform CAP analysis / analogue

As per the help files for ?capscale and the related ?dbrda, if we use Euclidean distances, as we do here, then capscale is identical to rda with the latter also being more efficient. Depending on the details of what one wants to achieve with this analysis which has its origins in PRIMER, you may want to explore more about these alternative functions.

# fit rda model from the vegan package
CAP_like_analysis <- vegan::rda( Y_muscle  ~ Species, 
                 data = mydata, 
                 dist = "euclidean")

# try to predict a value for an "Arrowtooth.Flounder"
# for example.
predict(CAP_like_analysis, 
        newdata = data.frame(Species = 
                               "Arrowtooth.Flounder"))

  Muscle.d13C Muscle.d15N
1   -20.07823    14.64197